Environmental data

0.1 Setup - load packages

library(ggplot2)
library(cowplot)
library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(colorBlindness)
library(scales)
library(bestNormalize)
library(Rmisc)

## Loading required package: lattice

## Loading required package: plyr

## ------------------------------------------------------------------------------

## You have loaded plyr after dplyr - this is likely to cause problems.
## If you need functions from both plyr and dplyr, please load plyr first, then dplyr:
## library(plyr); library(dplyr)

## ------------------------------------------------------------------------------

## 
## Attaching package: 'plyr'

## The following objects are masked from 'package:dplyr':
## 
##     arrange, count, desc, failwith, id, mutate, rename, summarise,
##     summarize

library(ggpubr)

## 
## Attaching package: 'ggpubr'

## The following object is masked from 'package:plyr':
## 
##     mutate

## The following object is masked from 'package:cowplot':
## 
##     get_legend

1 Substrate composition

#home comp
#sub.tah <- read.csv('/Volumes/Macintosh HD/Users/nicolakriefall/Google Drive (nicfall@bu.edu)/Moorea_revisions/env_revised/substrate.csv',header=TRUE)
#work comp
sub.tah <- read.csv('~/nicfall@bu.edu - Google Drive/My Drive/Moorea_revisions/env_revised/substrate.csv',header=TRUE)

sub <- subset(sub.tah,Marine.Area!="Tahiti Faaa")

sub$Substrate <- as.factor(sub$Substrate)
#str(sub)
##total proportion
sub$prop.test <- (sub$proportion)/3

##rename sites & habitats
sub$Marine.Area <- as.factor(sub$Marine.Area)
levels(sub$Marine.Area)

## [1] "Maatea"  "Tiahura"

levels(sub$Marine.Area) <- c("Mo'orea SE","Mo'orea NW")
levels(sub$Marine.Area)

## [1] "Mo'orea SE" "Mo'orea NW"

sub$Marine.Area <- factor(sub$Marine.Area, levels=c("Mo'orea NW","Mo'orea SE")) #reorder

sub$Habitat <- as.factor(sub$Habitat)
levels(sub$Habitat)

## [1] "Barrier reef" "Outer slope"

levels(sub$Habitat) <- c("Back","Fore")
levels(sub$Habitat)

## [1] "Back" "Fore"

##only interested in 2013 - year of collection
sub13 <- subset(sub,Year=="2013")
sub14 <- subset(sub,Year=="2014")

1.1 Plots

1.1.1 Figure - all the categories

ggplot(data=sub13, aes(x=Habitat, y=prop.test, fill=Substrate)) +
  geom_bar(stat="identity")+
  facet_wrap(~Marine.Area)+
  theme_cowplot()+
  ylab('Relative proportion')+
  xlab('Reef zone')

ggplot(data=sub14, aes(x=Habitat, y=prop.test, fill=Substrate)) +
  geom_bar(stat="identity")+
  facet_wrap(~Marine.Area)+
  theme_cowplot()+
  ylab('Relative proportion')+
  xlab('Reef zone')

#ggsave("substrate.png")

1.1.2 Figure - summed categories

sub13$Substrate <- gsub(" ","_",sub13$Substrate)

sub13$sub_cats <- recode(sub13$Substrate, Acropora="Coral", Asparagopsis="Algae", Cyanophyceae="Algae", Dead_coral="Abiotic substrata", Dictyota="Algae", Halimeda="Algae", Macroalgae="Algae", Millepora="Coral", Montastrea="Coral", Montipora="Coral", Napopora="Coral", Pavement="Abiotic substrata", Pavona="Coral", Pocillopora="Coral", Porites="Coral", Psammocora="Coral", Rubble="Abiotic substrata", Sand="Abiotic substrata", Stegastes_Turf="Algae", Synarea="Coral", Turbinaria="Coral")

pal <- hcl.colors(12,"Tofino")
show_col(pal)

colorz <- c("#3B415E","#8FC67C","#818CC0")

sub3 <- ggplot(data=sub13, aes(x=Habitat, y=prop.test, fill=sub_cats)) +
  geom_bar(stat="identity")+
  facet_wrap(~Marine.Area)+ #just need to remove 'nrow' here to go back to grid
  theme_cowplot()+
  ylab('Relative proportion')+
  xlab('Reef zone')+
  scale_fill_manual(values=colorz,name="Substrate",labels=c("Abiotic","Algae","Coral"))

sub3

#cvdPlot(sub3)

#ggsave("substrate_3cats.png",height=3.5) #can also go back to default if square

sub14$Substrate <- gsub(" ","_",sub14$Substrate)

sub14$sub_cats <- recode(sub14$Substrate, Acropora="Coral", Asparagopsis="Algae", Cyanophyceae="Algae", Dead_coral="Abiotic substrata", Dictyota="Algae", Halimeda="Algae", Macroalgae="Algae", Millepora="Coral", Montastrea="Coral", Montipora="Coral", Napopora="Coral", Pavement="Abiotic substrata", Pavona="Coral", Pocillopora="Coral", Porites="Coral", Psammocora="Coral", Rubble="Abiotic substrata", Sand="Abiotic substrata", Stegastes_Turf="Algae", Synarea="Coral", Turbinaria="Coral")

pal <- hcl.colors(12,"Tofino")
show_col(pal)

colorz <- c("#3B415E","#8FC67C","#818CC0")

sub4 <- ggplot(data=sub14, aes(x=Habitat, y=prop.test, fill=sub_cats)) +
  geom_bar(stat="identity")+
  facet_wrap(~Marine.Area)+ #just need to remove 'nrow' here to go back to grid
  theme_cowplot()+
  ylab('Relative proportion')+
  xlab('Reef zone')+
  scale_fill_manual(values=colorz,name="Substrate",labels=c("Abiotic","Algae","Coral"))

sub4

1.1.3 Figure - medium categories

sub13$sub_cats_med <- recode(sub13$Substrate, Acropora="Coral - Acropora", Asparagopsis="Algae - Asparagopsis", Cyanophyceae="Algae - Cyanophyceae", Dead_coral="Abiotic (e.g. rubble, sand)", Dictyota="Algae - Dictyota", Halimeda="Algae - Halimeda", Macroalgae="Algae - Macroalgae", Millepora="Coral - Millepora", Montastrea="Coral - Montastrea", Montipora="Coral - Montipora", Napopora="Coral - Napopora", Pavement="Abiotic (e.g. rubble, sand)", Pavona="Coral - Pavona", Pocillopora="Coral - Pocillopora", Porites="Coral - Porites", Psammocora="Coral - Psammocora", Rubble="Abiotic (e.g. rubble, sand)", Sand="Abiotic (e.g. rubble, sand)", Stegastes_Turf="Algae - Stegastes Turf", Synarea="Coral - Synarea", Turbinaria="Coral - Turbinaria")  

#greens for algae
#pal.gr <- hcl.colors(5,"Greens")
#show_col(pal.gr)
#pal.gr
#"#004616" "#30893B" "#81C07A" "#CAE8C1" "#F6FBF4"
#pal.bl <- hcl.colors(9,"BuPu")
#show_col(pal.bl)
#pal.bl
#"#540046" "#6E2072" "#8346A1" "#8674B8" "#9099CA" "#A3BAD9" "#BED6E6" "#DBECF3" "#F2FBFC"

colorz2 <- c("grey","#004616","#30893B","#81C07A","#CAE8C1","#F6FBF4","#540046","#6E2072","#8346A1","#8674B8","#9099CA","#A3BAD9","#BED6E6","#DBECF3","#F2FBFC")

sub.med <- ggplot(data=sub13, aes(x=Habitat, y=prop.test, fill=sub_cats_med)) +
  geom_bar(stat="identity")+
  facet_wrap(~Marine.Area)+ #just need to remove 'nrow' here to go back to grid
  theme_cowplot()+
  ylab('Relative proportion')+
  xlab('Reef zone')+
 scale_fill_manual(values=colorz2,name="Substrate")

sub.med

#ggsave("substrate_medcats.png",height=4) #can also go back to default if square
#cvdPlot(sub.med)

sub14$sub_cats_med <- recode(sub14$Substrate, Acropora="Coral - Acropora", Asparagopsis="Algae - Asparagopsis", Cyanophyceae="Algae - Cyanophyceae", Dead_coral="Abiotic (e.g. rubble, sand)", Dictyota="Algae - Dictyota", Halimeda="Algae - Halimeda", Macroalgae="Algae - Macroalgae", Millepora="Coral - Millepora", Montastrea="Coral - Montastrea", Montipora="Coral - Montipora", Napopora="Coral - Napopora", Pavement="Abiotic (e.g. rubble, sand)", Pavona="Coral - Pavona", Pocillopora="Coral - Pocillopora", Porites="Coral - Porites", Psammocora="Coral - Psammocora", Rubble="Abiotic (e.g. rubble, sand)", Sand="Abiotic (e.g. rubble, sand)", Stegastes_Turf="Algae - Stegastes Turf", Synarea="Coral - Synarea", Turbinaria="Coral - Turbinaria")  

#greens for algae
pal.gr <- hcl.colors(5,"Greens")
show_col(pal.gr)

pal.gr

## [1] "#004616" "#30893B" "#81C07A" "#CAE8C1" "#F6FBF4"

#"#004616" "#30893B" "#81C07A" "#CAE8C1" "#F6FBF4"
pal.bl <- hcl.colors(9,"BuPu")
show_col(pal.bl)

pal.bl

## [1] "#540046" "#6E2072" "#8346A1" "#8674B8" "#9099CA" "#A3BAD9" "#BED6E6"
## [8] "#DBECF3" "#F2FBFC"

#"#540046" "#6E2072" "#8346A1" "#8674B8" "#9099CA" "#A3BAD9" "#BED6E6" "#DBECF3" "#F2FBFC"

colorz2 <- c("grey","#004616","#30893B","#81C07A","#CAE8C1","#F6FBF4","#540046","#6E2072","#8346A1","#8674B8","#9099CA","#A3BAD9","#BED6E6","#DBECF3","#F2FBFC")

sub.med2 <- ggplot(data=sub14, aes(x=Habitat, y=prop.test, fill=sub_cats_med)) +
  geom_bar(stat="identity")+
  facet_wrap(~Marine.Area)+ #just need to remove 'nrow' here to go back to grid
  theme_cowplot()+
  ylab('Relative proportion')+
  xlab('Reef zone')+
 scale_fill_manual(values=colorz2,name="Substrate")

sub.med2

#ggsave("substrate_medcats.png",height=4) #can also go back to default if square
#cvdPlot(sub.med)

1.2 Synthesizing

ggarrange(sub.med,sub.med2,common.legend=TRUE,legend="right")

2 Nutrients

#home comp
#nuts <- read.csv("/Volumes/Macintosh HD/Users/nicolakriefall/Google Drive (nicfall@bu.edu)/Moorea_revisions/env_revised/nuts_13-14_renamed.csv",header=T)
#work comp
nuts <- read.csv("~/nicfall@bu.edu - Google Drive/My Drive/Moorea_revisions/env_revised/nuts_13-14_renamed.csv",header=T)

str(nuts)

## 'data.frame':    48 obs. of  11 variables:
##  $ year             : int  2013 2013 2013 2013 2013 2013 2013 2013 2013 2013 ...
##  $ month            : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ date             : chr  "1/20/13" "2/20/13" "3/20/13" "4/20/13" ...
##  $ site             : chr  "Tiahura" "Tiahura" "Tiahura" "Tiahura" ...
##  $ habitat          : chr  "Fore reef" "Fore reef" "Fore reef" "Fore reef" ...
##  $ replicate        : int  3 3 3 3 3 3 3 3 3 3 ...
##  $ Phosphate_P2O5_uM: num  0.527 0.47 0.258 0.293 0.468 0.156 0.031 0.175 0.142 0.319 ...
##  $ Nitrates_NO3_uM  : num  0.108 0.101 0.114 0.164 0.023 0.186 0.039 0.049 0.009 0.011 ...
##  $ Nitrites_NO2_uM  : num  0.069 0.008 0.028 0.036 0.013 0.006 0.019 0.002 0.034 0.005 ...
##  $ Silice_SiO2_uM   : num  2.58 7.11 1.28 1.08 1.2 ...
##  $ Ammonium_NH4_uM  : num  4.637 2.198 1.06 0.473 0.55 ...

2.1 From revisions - by month

2.1.1 Phosphate revised

nuts$ym <- paste0(nuts$year,"_",nuts$month)
nuts$date2 <- as.Date(nuts$date,"%m/%d/%y")
str(nuts)

## 'data.frame':    48 obs. of  13 variables:
##  $ year             : int  2013 2013 2013 2013 2013 2013 2013 2013 2013 2013 ...
##  $ month            : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ date             : chr  "1/20/13" "2/20/13" "3/20/13" "4/20/13" ...
##  $ site             : chr  "Tiahura" "Tiahura" "Tiahura" "Tiahura" ...
##  $ habitat          : chr  "Fore reef" "Fore reef" "Fore reef" "Fore reef" ...
##  $ replicate        : int  3 3 3 3 3 3 3 3 3 3 ...
##  $ Phosphate_P2O5_uM: num  0.527 0.47 0.258 0.293 0.468 0.156 0.031 0.175 0.142 0.319 ...
##  $ Nitrates_NO3_uM  : num  0.108 0.101 0.114 0.164 0.023 0.186 0.039 0.049 0.009 0.011 ...
##  $ Nitrites_NO2_uM  : num  0.069 0.008 0.028 0.036 0.013 0.006 0.019 0.002 0.034 0.005 ...
##  $ Silice_SiO2_uM   : num  2.58 7.11 1.28 1.08 1.2 ...
##  $ Ammonium_NH4_uM  : num  4.637 2.198 1.06 0.473 0.55 ...
##  $ ym               : chr  "2013_1" "2013_2" "2013_3" "2013_4" ...
##  $ date2            : Date, format: "2013-01-20" "2013-02-20" ...

gg.new.pho <- ggplot(nuts,aes(x=date2,y=Phosphate_P2O5_uM,color=habitat,group=habitat,shape=habitat))+
  geom_line()+
  theme_cowplot()+
  geom_point()+
  geom_vline(xintercept=as.numeric(as.Date("2013-07-26")),linetype=4)+
  geom_vline(xintercept=as.numeric(as.Date("2013-08-09")),linetype=4)+
  scale_shape_manual(values=c(16,15))+
  scale_colour_manual(values=c("#ED7953FF","#8405A7FF"),labels=c("Back reef","Fore reef"))+
  #theme(legend.position="none")+
  xlab("")+
  ylab("Phosphates (µM)")+
  labs(color="Habitat",shape="Habitat")
gg.new.pho

2.1.2 Nitrates revised

gg.new.nia <- ggplot(nuts,aes(x=date2,y=Nitrates_NO3_uM,color=habitat,group=habitat,shape=habitat))+
  geom_line()+
  theme_cowplot()+
  geom_point()+
  geom_vline(xintercept=as.numeric(as.Date("2013-07-26")),linetype=4)+
  geom_vline(xintercept=as.numeric(as.Date("2013-08-09")),linetype=4)+
  scale_shape_manual(values=c(16,15))+
  scale_colour_manual(values=c("#ED7953FF","#8405A7FF"),labels=c("Back reef","Fore reef"))+
  #theme(legend.position="none")+
  xlab("")+
  ylab("Nitrates (µM)")+
  labs(color="Habitat",shape="Habitat")
gg.new.nia

2.1.3 Nitrites revised

gg.new.nii <- ggplot(nuts,aes(x=date2,y=Nitrites_NO2_uM,color=habitat,group=habitat,shape=habitat))+
  geom_line()+
  theme_cowplot()+
  geom_point()+
  geom_vline(xintercept=as.numeric(as.Date("2013-07-26")),linetype=4)+
  geom_vline(xintercept=as.numeric(as.Date("2013-08-09")),linetype=4)+
  scale_shape_manual(values=c(16,15))+
  scale_colour_manual(values=c("#ED7953FF","#8405A7FF"),labels=c("Back reef","Fore reef"))+
  #theme(legend.position="none")+
  xlab("")+
  ylab("Nitrites (µM)")+
  labs(color="Habitat",shape="Habitat")
gg.new.nii

2.1.4 Silica revised

gg.new.sil <- ggplot(nuts,aes(x=date2,y=Silice_SiO2_uM,color=habitat,group=habitat,shape=habitat))+
  geom_line()+
  theme_cowplot()+
  geom_point()+
  geom_vline(xintercept=as.numeric(as.Date("2013-07-26")),linetype=4)+
  geom_vline(xintercept=as.numeric(as.Date("2013-08-09")),linetype=4)+
  scale_shape_manual(values=c(16,15))+
  scale_colour_manual(values=c("#ED7953FF","#8405A7FF"),labels=c("Back reef","Fore reef"))+
  #theme(legend.position="none")+
  xlab("")+
  ylab("Silica (µM)")+
  labs(color="Habitat",shape="Habitat")
gg.new.sil

2.1.5 Ammonium revised

gg.new.amm <- ggplot(nuts,aes(x=date2,y=Ammonium_NH4_uM,color=habitat,group=habitat,shape=habitat))+
  geom_line()+
  theme_cowplot()+
  geom_point()+
  geom_vline(xintercept=as.numeric(as.Date("2013-07-26")),linetype=4)+
  geom_vline(xintercept=as.numeric(as.Date("2013-08-09")),linetype=4)+
  scale_shape_manual(values=c(16,15))+
  scale_colour_manual(values=c("#ED7953FF","#8405A7FF"),labels=c("Back reef","Fore reef"))+
  #theme(legend.position="none")+
  xlab("")+
  ylab("Ammonium (µM)")+
  labs(color="Habitat",shape="Habitat")
gg.new.amm

2.2 Results summarized - revised

ggarrange(gg.new.pho,gg.new.nia,gg.new.nii,gg.new.sil,gg.new.amm,labels="AUTO",common.legend=TRUE,legend="right",nrow=5)

#ggsave("env_nutrients_revised.pdf",height=12)

2.3 Box plots

2.3.1 Wrap by year

ggplot(nuts,aes(x=habitat,y=Phosphate_P2O5_uM))+
  geom_boxplot()+
  facet_wrap(~year)

ggplot(nuts,aes(x=habitat,y=Nitrates_NO3_uM))+
  geom_boxplot()+
  facet_wrap(~year)

ggplot(nuts,aes(x=habitat,y=Nitrites_NO2_uM))+
  geom_boxplot()+
  facet_wrap(~year)

ggplot(nuts,aes(x=habitat,y=Silice_SiO2_uM))+
  geom_boxplot()+
  facet_wrap(~year)

ggplot(nuts,aes(x=habitat,y=Ammonium_NH4_uM))+
  geom_boxplot()+
  facet_wrap(~year)

2.4 Violin plots

2.4.1 Phosphates

gg.pho <- ggplot(nuts,aes(x=habitat,y=Phosphate_P2O5_uM,color=habitat,shape=habitat))+
  geom_violin()+
  #geom_jitter(alpha=0.7)+
  theme_cowplot()+
  geom_boxplot(width=0.1,alpha=0)+
  scale_shape_manual(values=c(16,15),labels=c("BR","FR"))+
  scale_x_discrete(labels=c("BR","FR"))+
  scale_colour_manual(values=c("#ED7953FF","#8405A7FF"),labels=c("Back reef","Fore reef"))+
  theme(legend.position="none")+
  xlab("Habitat")+
  ylab("Phosphates (µM)")
gg.pho

library(car)

## Loading required package: carData

## 
## Attaching package: 'car'

## The following object is masked from 'package:dplyr':
## 
##     recode

nuts$habitat <- as.factor(nuts$habitat)
shapiro.test(log(nuts$Phosphate_P2O5_uM))

## 
##  Shapiro-Wilk normality test
## 
## data:  log(nuts$Phosphate_P2O5_uM)
## W = 0.95694, p-value = 0.07596

leveneTest(nuts$Phosphate_P2O5_uM,nuts$habitat)#ns

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  1  0.0922 0.7627
##       46

summary(aov(Phosphate_P2O5_uM~habitat,data=nuts)) #ns

##             Df Sum Sq Mean Sq F value Pr(>F)
## habitat      1 0.0025 0.00251   0.049  0.826
## Residuals   46 2.3566 0.05123

#wilcox.test(Phosphate_P2O5_uM~habitat,data=nuts,exact=FALSE) #ns

2.4.2 Nitrates

gg.nia <- ggplot(nuts,aes(x=habitat,y=Nitrates_NO3_uM,color=habitat,shape=habitat))+
  geom_violin()+
  #geom_jitter(alpha=0.7)+
  geom_boxplot(width=0.1,alpha=0)+
  theme_cowplot()+
  scale_shape_manual(values=c(16,15),labels=c("BR","FR"))+
  scale_x_discrete(labels=c("BR","FR"))+
  scale_colour_manual(values=c("#ED7953FF","#8405A7FF"),labels=c("Back reef","Fore reef"))+
  theme(legend.position="none")+
  xlab("Habitat")+
  ylab("Nitrates (µM)")
gg.nia

shapiro.test(log(nuts$Nitrates_NO3_uM)) #no

## 
##  Shapiro-Wilk normality test
## 
## data:  log(nuts$Nitrates_NO3_uM)
## W = 0.91919, p-value = 0.002784

leveneTest(nuts$Nitrates_NO3_uM,nuts$habitat)#ns

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  1  2.3953 0.1286
##       46

#summary(aov(Nitrates_NO3_uM~habitat,data=nuts)) #p < 0.01**
wilcox.test(Nitrates_NO3_uM~habitat,data=nuts,exact=FALSE) #p < 0.01**

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  Nitrates_NO3_uM by habitat
## W = 435, p-value = 0.002519
## alternative hypothesis: true location shift is not equal to 0

# Wilcoxon rank sum test with continuity correction
# 
# data:  Nitrates_NO3_uM by habitat
# W = 435, p-value = 0.002519

sum.nia <- summarySE(data=nuts,measurevar="Nitrates_NO3_uM",groupvars="habitat")
sum.nia

##     habitat  N Nitrates_NO3_uM         sd         se         ci
## 1 Back reef 24        0.171250 0.08852278 0.01806964 0.03737989
## 2 Fore reef 24        0.094625 0.05985476 0.01221780 0.02527445

2.4.3 Nitrites

gg.nii <- ggplot(nuts,aes(x=habitat,y=Nitrites_NO2_uM,color=habitat,shape=habitat))+
  geom_violin()+
  #geom_jitter(alpha=0.7)+
  geom_boxplot(width=0.1,alpha=0)+
  theme_cowplot()+
  scale_shape_manual(values=c(16,15),labels=c("BR","FR"))+
  scale_x_discrete(labels=c("BR","FR"))+
  scale_colour_manual(values=c("#ED7953FF","#8405A7FF"),labels=c("Back reef","Fore reef"))+
  theme(legend.position="none")+
  xlab("Habitat")+
  ylab("Nitrites (µM)")
gg.nii

shapiro.test(nuts$Nitrites_NO2_uM)

## 
##  Shapiro-Wilk normality test
## 
## data:  nuts$Nitrites_NO2_uM
## W = 0.94972, p-value = 0.03897

#summary(aov(Nitrites_NO2_uM~habitat,data=nuts)) #p < 0.05*
wilcox.test(Nitrites_NO2_uM~habitat,data=nuts,exact=FALSE) #p < 0.01**

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  Nitrites_NO2_uM by habitat
## W = 417, p-value = 0.008029
## alternative hypothesis: true location shift is not equal to 0

sum.nii <- summarySE(data=nuts,measurevar="Nitrites_NO2_uM",groupvars="habitat")
sum.nii

##     habitat  N Nitrites_NO2_uM         sd          se          ci
## 1 Back reef 24      0.03216667 0.01396476 0.002850544 0.005896801
## 2 Fore reef 24      0.02083333 0.01636185 0.003339849 0.006909003

2.4.4 Silica

gg.sil <- ggplot(nuts,aes(x=habitat,y=Silice_SiO2_uM,color=habitat,shape=habitat))+
  geom_violin()+
  #geom_jitter(alpha=0.5)+
  geom_boxplot(width=0.1,alpha=0)+
  theme_cowplot()+
  scale_shape_manual(values=c(16,15),labels=c("BR","FR"))+
  scale_x_discrete(labels=c("BR","FR"))+
  scale_colour_manual(values=c("#ED7953FF","#8405A7FF"),labels=c("Back reef","Fore reef"))+
  theme(legend.position="none")+
  xlab("Habitat")+
  ylab("Silica (µM)")
gg.sil

summary(aov(Silice_SiO2_uM~habitat,data=nuts)) #ns

##             Df Sum Sq Mean Sq F value Pr(>F)
## habitat      1   0.19  0.1948   0.161   0.69
## Residuals   46  55.59  1.2084

wilcox.test(Silice_SiO2_uM~habitat,data=nuts) #ns

## 
##  Wilcoxon rank sum exact test
## 
## data:  Silice_SiO2_uM by habitat
## W = 291, p-value = 0.9593
## alternative hypothesis: true location shift is not equal to 0

2.4.5 Ammonium

gg.amm <- ggplot(nuts,aes(x=habitat,y=Ammonium_NH4_uM,color=habitat,shape=habitat))+
  geom_violin()+
  #geom_jitter(alpha=0.5)+
  geom_boxplot(width=0.1,alpha=0)+
  theme_cowplot()+
  scale_shape_manual(values=c(16,15),labels=c("BR","FR"))+
  scale_x_discrete(labels=c("BR","FR"))+
  scale_colour_manual(values=c("#ED7953FF","#8405A7FF"),labels=c("Back reef","Fore reef"))+
  theme(legend.position="none")+
  xlab("Habitat")+
  ylab("Ammonium (µM)")
gg.amm

summary(aov(Ammonium_NH4_uM~habitat,data=nuts)) #ns

##             Df Sum Sq Mean Sq F value Pr(>F)
## habitat      1   0.00  0.0012   0.001  0.974
## Residuals   46  50.02  1.0873

wilcox.test(Ammonium_NH4_uM~habitat,data=nuts) #ns

## 
##  Wilcoxon rank sum exact test
## 
## data:  Ammonium_NH4_uM by habitat
## W = 333, p-value = 0.3622
## alternative hypothesis: true location shift is not equal to 0

2.5 Multipanel violin plot

Results: - phosphates not sig, log-transformed anova - nitrates - nitrites - silica - ammonium

library(ggpubr)

ggarrange(gg.pho,gg.nia,gg.nii,gg.sil,gg.amm,labels="AUTO")

#ggsave("env_nutrients.pdf",width=7)

3 Flow

#home comp
#flowdat <- read.csv("/Volumes/Macintosh HD/Users/nicolakriefall/Google Drive (nicfall@bu.edu)/Moorea_revisions/env_revised/flow_rates_reformatted.csv",header=T)
#work comp
flowdat <- read.csv("~/nicfall@bu.edu - Google Drive/My Drive/Moorea_revisions/env_revised/flow_rates_reformatted.csv",header=T)

ggplot(flowdat,aes(x=site,y=speed_ms,fill=site))+
  geom_boxplot()+
  theme_cowplot()+
  scale_shape_manual(values=c(16,15))+
  scale_colour_manual(values=c("#ED7953FF","#8405A7FF"))+
  scale_fill_manual(values=c("#ED7953FF","#8405A7FF"))+
  ylab("Flow (m/s)")+
  xlab("Site-reef zone")+
  theme(legend.position="none")

#ggsave("mnw.flow.pdf",width=2.5,height=2)

best.speed <- bestNormalize(flowdat$speed_ms)

flowdat$speed.t <- best.speed$x.t

a.flow <- aov(speed.t~site,data=flowdat)
summary(a.flow)

##                Df Sum Sq Mean Sq F value Pr(>F)    
## site            1   6128    6128   11195 <2e-16 ***
## Residuals   13536   7409       1                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

plot(a.flow)

## hat values (leverages) are all = 0.0001477323
##  and there are no factor predictors; no plot no. 5

flow.sd <- summarySE(data=flowdat,measurevar="speed_ms",groupvars="site")
flow.sd

##    site    N   speed_ms         sd           se           ci
## 1 MNW-B 6769 0.06481833 0.02387819 0.0002902279 0.0005689381
## 2 MNW-F 6769 0.10730170 0.02182916 0.0002653229 0.0005201163

sessionInfo()

## R version 4.1.0 (2021-05-18)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Mojave 10.14.6
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] car_3.0-12           carData_3.0-4        ggpubr_0.4.0        
##  [4] Rmisc_1.5            plyr_1.8.6           lattice_0.20-45     
##  [7] bestNormalize_1.8.2  scales_1.1.1         colorBlindness_0.1.9
## [10] dplyr_1.0.7          cowplot_1.1.1        ggplot2_3.3.5       
## 
## loaded via a namespace (and not attached):
##  [1] tidyr_1.1.4        splines_4.1.0      foreach_1.5.1      prodlim_2019.11.13
##  [5] assertthat_0.2.1   highr_0.9          doRNG_1.8.2        yaml_2.2.1        
##  [9] globals_0.14.0     ipred_0.9-12       pillar_1.6.4       backports_1.3.0   
## [13] glue_1.5.0         digest_0.6.28      ggsignif_0.6.3     colorspace_2.0-2  
## [17] recipes_0.1.17     htmltools_0.5.2    Matrix_1.3-4       timeDate_3043.102 
## [21] pkgconfig_2.0.3    broom_0.7.10       listenv_0.8.0      purrr_0.3.4       
## [25] gower_0.2.2        lava_1.6.10        tibble_3.1.6       farver_2.1.0      
## [29] generics_0.1.1     usethis_2.1.3      ellipsis_0.3.2     withr_2.4.2       
## [33] nnet_7.3-16        survival_3.2-13    magrittr_2.0.1     crayon_1.4.2      
## [37] evaluate_0.14      fs_1.5.0           fansi_0.5.0        future_1.23.0     
## [41] parallelly_1.28.1  doParallel_1.0.16  MASS_7.3-54        rstatix_0.7.0     
## [45] class_7.3-19       tools_4.1.0        lifecycle_1.0.1    stringr_1.4.0     
## [49] munsell_0.5.0      rngtools_1.5.2     compiler_4.1.0     jquerylib_0.1.4   
## [53] gridGraphics_0.5-1 rlang_0.4.12       grid_4.1.0         iterators_1.0.13  
## [57] labeling_0.4.2     rmarkdown_2.11     gtable_0.3.0       codetools_0.2-18  
## [61] abind_1.4-5        DBI_1.1.1          R6_2.5.1           gridExtra_2.3     
## [65] lubridate_1.8.0    knitr_1.36         fastmap_1.1.0      future.apply_1.8.1
## [69] utf8_1.2.2         nortest_1.0-4      butcher_0.1.5      stringi_1.7.5     
## [73] parallel_4.1.0     Rcpp_1.0.7         vctrs_0.3.8        rpart_4.1-15      
## [77] tidyselect_1.1.1   xfun_0.28